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The numerical results clearly show how the Stokes flow behavior for low Reynolds 
numbers, and the boundary layer behavior for high Reynolds numbers, are approached in 
the appropriate limits. Besides showing the flow streamlines, results have been presented 
for the torque and the skin friction behavior. It is shown that the present results are in 
excellent agreement with both available experimental data, and previously obtained 
numerical data. The radial equatorial jet which develops with increasing Reynolds numbers 
has been observed as expected from boundary layer collision behavior. No separation was 
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I. INTRODUCTION 


The study of direct numerical simulation of steady viscous flow due to a sphere 
which rotates about a diameter with constant angular velocity in a fluid at rest has 
historically been the focus of a great deal of scientific attention. Most of the previous 
studies were based on finite-difference schemes. 

In this study, a different numerical scheme based on spectral methods was used. 
The validity of the numerical method and the results were checked with the previous 
studies. Spectral methods have become increasingly popular in recent years, especially 
since the development of fast transform methods, with applications in numerical weather 
prediction, numerical simulation of laminar and turbulent flows and other problems where 
high accuracy is desired for complicated results. 

The objective of this study is to compute the numerical simulation of a steady flow 
about a spinning sphere by using the spectral numerical methods. A Matlab program was 
developed by using spectral methods. However the method itself was not compared to the 
other methods used in previous studies. The reliability and validity of the numerical 
method were investigated to improve our current ability to compute and predict the 


evolution of flow for a particular geometry and conditions. 





Il. BACKGROUND STUDIES 


The spectral method is based on the collocation method which appears to have 
been used first by Slater and Kantorovic in specific applications in 1934. It was developed 
as a general method for solving ordinary differential equations. This method was revived 
by Wrigt in 1964. The applications of Chebyshev polynomial expansions to the initial 
value problems were involved in these studies. 

The spectral collocation method was applied to partial differential equations for 
spatially periodic problems by Kreiss and Oliger, who called it the Fourier method, and 
Orszag, who termed it “pseudospectral’”. These were the earliest applications of spectral 
collocation or the pseudospectral method. This approach was very attractive because of 
its application to variable-coefficient and even non-linear problems. 

The Galerkin approach which depends on the same trial and test functions, was 
applied to a meteorological model by Silberman in 1954. This was the first serious 
application of spectral methods to PDEs. In 1970, Machenhauer and Rasmussen 
developed transform methods for evaluating convolution sums arising from quadratic non- 
linearity. Spectral Galerkin methods became practical for high resolution calculations of 
such non-linear problems. Applications in fluid dynamics were reviewed in the symposium 
proceedings edited by Voigt, Gottlieb and Hussaini 1n1984. 

For the problem of the flow induced by a spinning sphere being considered in this 
study, the flow behavior for small values of the Reynolds number, Stokes [Ref. 1], Lamb 
[Ref. 2], Bickley [Ref .3], Collins [Ref. 4], Thomas and Walters [Ref. 5], Ovseenko 
[Ref. 6], and Takagi [Ref. 7] have given theoretical results. For high Reynolds number, the 
boundary-layer equations have been studied by Howarth [Ref. 8], Nigam [Ref. 9], 
Stewartson [Ref. 10], Fox [Ref. 11], Banks [Refs. 12 & 13], Manohar [Ref. 14], and 
Singh [Ref. 15]. Over a large range of values of Reynolds number the flow has been 
investigated experimentally by Kobashi [Ref. 16], Bowden and Lord [Ref. 17], Kreith et 


al [Ref. 18], and Sawatzki [Ref. 19]. Another numerical method based on the use of 


specialized techniques to obtain an approximation of the finite-difference form of the full 
partial differential equations was performed by Allen and Southwell [Ref. 20], Dennis 
[Refs. 21 & 22], Roscoe [Refs. 23 & 24], Spalding [Ref. 25] and Dennis, Ingham and 
Singh [Ref. 26]. 


Ill. GOVERNING EQUATIONS 


A. DERIVATION OF EQUATIONS 
The governing equations for the problem are as follow; 


au 





, 
> +(u.V)u= a Vp+vV-u (3.1) Momentum equation 
V.u=0 (3.2) Continuity equation 
Define dimensionless variables; 
r . wu . Vp a 
6 a © U ; 0U? a 


where a : radius of the sphere, 
U : free stream velocity, 
Oo : density of the fluid, 
Vv : kinematic viscosity of the fluid. 


Introduce the dimensionless variables into Equation (3.1); 
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Uv ou UU? . 1U" U 
cade cet "Vv i oe V * —\j* * ; 
a ae _ )u Da Das eo u (3.3) 
Divide Equation (3.3) by a 
nL 
ou- alley yu as ‘wy (3.4) 
—+—(u Vj =-— : 
ot Vv V e 


Hereafter the superscript * will be omitted for nondimensional terms. Define the 


U(2a) 


Reynolds number based on the diameter, Reg = , and take the Curl of Equation 


(3.4) to be able to eliminate the pressure term and introduce vorticity, @ = V xu_ where 
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(WO=1,0),+1eWo+lo Wo, 


The transport equation for the ®-component of the vorticity is 


6) , 


; : 2 i : 
a Coie ~VXUXMsie) = Re vy “Tat Mole (3.5) 


Convert the velocity in Equation (3.5) by defining the axisymmetric Stream function, y, 


and an “angular circulation”, Q2. 
V —_ a Azimuthal direct 3.6 
= X\(orO.OC-— + — : 
u <A). | ae ® ; Azimuthal direction (3.6) 
In spherical coordinates, Equation (3.6) is obtained as 


(| 8), |) 24 
"— r- sin@ 00 it rsin® or '@ ‘eho e oe 


ie ele tot Ua (328) 








where 


LS eee ee 
In the free stream y= 5 Ur sin’ 8 or y= ae sin” @ in nondimensional form. 


By using the definition of the stream function it can be shown that 


-2¥,t (2 wv )(—) 
a=] a r° sin8 S(4 fale) rsin®@ (3.9) 








Define operator D* as 











p= 2, ee <) (3.10) 
ar? rr? 0B sin 8 00 
to obtain the azimuthal component of the vorticity; 
n2 
= D’y (Sug 


Call [VX (u X @eie)].i@ in Equation (3.5) nonlinear Term G(r,0,t). By carrying out the 


calculations in spherical coordinates 


1 |ocy,D’yw) — 
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is obtained, where 
uu = cos® 
 wulare 
P (3.13) 
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Combine all terms to get 
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Equation (3.14) is called the Vorticity-Stream function form of the Momentum 


Equation in axisymmetric spherical coordinate system. For computational purposes, 


which will become clear later, define a modified stream function C(r,8,t), which is related 


to the usual Stokes stream function (y) by the relation; 


yv =rCsin®@ 


Also for any function f(r,®), define operator D~ such that 


D°(fr sin®@) = rsinO9D’f 











where 
] 
D* =V* -————_. 
r~ sin’ 9 
This gives 
@, =—-D°C 
= Vai.) + 
o (C1,) Aarons 
1 oa (Csin8) 
= —(Cs 
vl rsin9 00 ss 
12g 
where Ug = ie r 
Q 
u — 
° rsin® 


(Sais) 


(3.16) 


(Salgy 


3-13} 


The modified stream function C is written as the sum C aor ae where C=0. 


(which is the prescribed quiescent free stream condition), and c represents the disturbance 


produced by the presence of the spinning sphere. 
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| 7 Re 
Returning to Equation (3.5), rearrange by multiplying both sides with r° = 


Re eee 








Se 2n2 
> Ot r 5 G(r,0,t)+r°D°@ (3.19) 
r’D*C=-r’o (3.20) 
where 
G(r,0,t)=[Vx(ux@gig)| ig (3.21) 


Note that the operator D? is defined as 


ral 2) <(s o< | = 522) 
TWP esa a: r°> sin’ @ S 


Divide r°D* in two parts which are D? and D? respectively and obtain 


2 zy o 2 
PD =) 5b, Bevo 
sine Cae a 
where D? and Djare as defined in Chapter IV. 


The nonlinear term in Equation (3.21), G(r,6,t), consists of 8 terms. These are 
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Time integration of vorticity equation (3.19) is accomplished through the use of an 


explicit second-order Adams-Bashforth scheme for the nonlinear terms and an implicit 
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second-order Crank-Nicolson scheme for the viscous and linear terms. The calculations 
are made in physical] r-space and spectral 6-space. If we denote vorticity @(r,t) as the 


vector of N sine series then the discretized form of Equation (3.19) is 
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, Re Disar 7, > Re 
“| 7 Jer - (3.24) 
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where & and B are suitably chosen weighting parameters in the implicit scheme. 


Rearrange Equation (3.24) and obtain 








2 2 Rey _ 2 2 Reg 2 
1D are Ie iN Oe eer at ic O21 Re ,|aG, +BG,_,,| (3.23) 


where a= = p=-~ and G=(Vx (ux Wgig))-ig 


Ps + Alc(r,t + At) =—r?@(r,t + At) (3.26) 


We find the governing equation for ug 1n same way: 











dQ 61] A(¥,Q) a 
—+—| —-——_ |= DQ 3.27 
ot 2) ae Re, Ca 
h = d C 
W ——_— = = 
Oe wae rSin8 


and introduce D Q=D™ (u,rSin@) = rSin@D*u, and divide each side by (r sin8) 


If Equation (3.27) is rearranged, the result is 
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a 7 ao 5(r.0) Re 


The nonlinear term in Equation (3.28) is called H(r,9,t) and consists of 
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Ug iS written as the sum of two parts: potential and disturbance parts, 


= ! —  Sin@ 
Uy =Ue FUy where Ug =—>— Is the Stokes flow 
r 





solution which satisfies Dus = Q. The final form of Equation (3.28) is 


OU g 
ot 








-_?(u;) (3.29) 
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3 = (e+ a0) SE se sD pI " 20) ee 2) — [3H(t)-H(t—At)] (3.31) 
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B. BOUNDARY CONDITIONS 


The boundary conditions for the numerical scheme are 


= dc aC | 
c—==C=0 = 2S = ao. — 0 eer (3:33) 
or or 
c=O W=0 w=0 atr=r.. (3.34) 


The Neumann Boundary Conditions in Equation (3.33) are handled with the 
Green’s Function Method that will be explained later. Zero boundary conditions in the 6 


direction are automatically satisfied with the choice of the sine expansions. 


IV. METHOD OF SOLUTION 


A. NUMERICAL METHOD 

A pseudospectral or collocation method based on the technique developed by 
Marcus & Tuckerman [Ref. 29], is used to solve the equations of motion in time. Spectral 
methods may be viewed as an extreme development of weighted residual (MWR) which is 
generally used in a discretization scheme for differential methods. The trial functions and 
test functions are chosen as infinitely differentiable global functions. This is one of the 
features which distinguish spectral methods from finite-element and finite-difference 
methods. A pseudospectral transform method is used in this study. The approach taken in 
the transform method is to use the Fast Fourier Transform (FFT) to transform the 
functions to spectral space or Inverse Fast Fourier Transform (IFFT) to transform the 
equations to physical space. 

A physical process can be described either in time domain, by the values of the 
same h as a function of t, e.g., h(t), or else in frequency domain, where the process is 
specified by giving its amplitude, H, as a function of frequency f, that is H(f), with -00 <f < 
co, So once it is known whether the function is either odd or even, this information can be 
used in transformations. If h(t) 1s real and even then H(f) is real and even, if h(t) is real 
and odd then H(f) is imaginary and odd. It is necessary to go back and forth between 
these two domains by means of Fourier transform equations. As explained below, odd 
functions, using sine series expansions for @, c and ue, are used. The same properties are 
valid for derivatives of these expansions. Using these properties will increase 
computational efficiency. 

Functions are presented both in spectral space as a finite series of basis functions 
and by values at collocation grid points in physical space. The vorticity, w, the modified 
stream function, C, and the velocity in the ® direction, ug, are expanded as Chebyshev 
polynomials in the radial direction and as sine series in the transverse direction. The 


products and the derivatives in radia] direction are done in physical space while the 


I 


derivatives in transverse direction are obtained in spectral space. Since each sine term in 
the expansion satisfies the homogenous 9 boundary conditions (w = 0, ‘¥ = 0, ue =0 at 
0=0,7) exactly, no further 9 boundary conditions need to be applied. 


In general let a variable, say @, be expressed as 


co(r,0,t) = > @, (z, t)Sin(0, ) (4.1a) 
with 
@,(z,t)= 90, (t)T., (2) (4.1b) 
m=0 


where T,,(z) is Chebyshev Polynomial with -1 <z < 1 





a — a SO ON 4.2 

seis n= 1,2,... (4.2) 
mit 

2=Cof ™) m = OF (4.3) 
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An algebraic map is used to map the radial interval | <r <r.. to the interval 


-1 <z<1 where 








I+Z 
P= eb rellizlesel (4.4) 
b-z 
Ze - 
with b=1+ (Note that rp = 1, which 1s the surface of the sphere.) 
es 0 


L is a scaling parameter which is used to map the radial interval to the z-interval. 
A finite, but large outer r., was chosen to avoid regularity problems in the radial 
differentiation. 

The sine expansions (4.1a) for w, C and ug satisfy the homogenous @ boundary 
conditions and match exactly the periodicity and symmetry conditions in 8. The solution 
of the elliptic equation (3.20) is required to obtain the modified stream function, C, from 
vorticity @ at any time level. Following [Ref. 29], separable derivative operators, D, and 


2, ° 
De’, are defined as 
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r°-D* =D? +— 
sin’ 8 


aed 
ror or 


» fn Af. d 
D; =| sinc 56 sino nl 


In spectral 6-space the effect of operator Dg’ at any fixed radial location (physical r-space) 


D; (4.5) 
where 


D 
(4.6) 


can be written as 


2 





Dé. f1@)=|sin* 0-2 a 


) 
+ sin 9.cos 9— -— JE sin( 9) 


Mz 


=f, 145 P)sinGjoy +4 je sin(26 + j@) += jj sin(20-~ 8) (4.7) 


J 
lt 


If this expression is written in the form of D3f(0) =|s, |r, 


Gig) ial 
i a 1=jt+2 
—|-— i=] 
S, = 2 (4.8) 
(i+2)04+1) +5 
4 
0 otherwise 


is obtained. 


This matrix is NxN if the representation for D; f is truncated to N Fourier terms. 


Similarly, the operation of multiplying a function f{(@) by sin™(@) produces another 


NxN matrix which is called A matrix where the only non-zero terms are 


al(i) = —-i(i+1) = j 


i= 4. 
As eevee 1< j,1+)—cven ase 


is: 


Time integration of the equations are obtained by using the scheme which is 
explained in Chapter III. 

The nonlinear terms, G and H, have to be carefully handled by transforming the 
spectral coefficients to physical space first, then multiplying the terms in physical space, 
and finally transforming the result of these products back to spectral space. 

The radial derivatives are evaluated by collocation methods and expressed as 
matrix operations on the vector of function values at the grid points in physical r-space. 


For (M+1) collocation points, we introduce the (M+1)x(M+1) matrix operation, that is 


Dea) = 1B) +} acy-A2 0} (4.10) 


where, I is the identity matrix. The operator (D,-r°(Re/At)+A] in (3.25) and (3.32) may 


now be written in block matrix form as 


Ba 0 a2(1)I 0 
0 Dee) 0 a2(2)I 
0 0 Dee) 0 
0 0 0 D‘” (4) 


(4.11) 


This matrix acts on the vector of a length of (M+1I)xN for vorticity coefficients. 
Equations (3.25) and (3.26) are upper triangular, block matrix problems in physical r- 
space and spectral @-space that can be solved by any solver with the Dirichlet and the 
Neumann boundary conditions discussed below. The equation (3.32) is treated in the 
same manner and leads to an upper triangular, block matrix problem. 

B. GREEN’S FUNCTION METHOD 
To enforce the radial boundary conditions in c-@ equations, (3.25) and (3.26), 


Green’s Function Method has been developed which gives the correct boundary 


16 


conditions, (3.33) and (3.34), and allows the surface vorticity to develop naturally. If @ 


and c consist of two parts, homogenous and particular solutions, the following is written: 


N 


N 
O=0, + 22,0, and C=C, + DAE j= 2a (4.12) 
j= F 


If (4.12) is introduced into (3.25) and (3.26), the homogenous parts as follow with 


corresponding boundary conditions occur. 


E*6 (1,8) = 0 (4.13) H*¢. = OF (4.14) 
@ (r= 1,0;) =9, c,(r =1,0,)=0 
(1 = ,0;) =0 C (r=0,0;)=0 


The particular parts with corresponding boundary conditions are 


E*@,(r,0,t + At)=R(r,6,t) (4.15) Hc, (r,0,t + At)=-r°o, (1,8, t) (4.16) 
@,(r =1,8) =0 c,(r=1,8) =—C.,_, 
QO, (r =o,8)=0 c, (r= 0,0) =0 


E’ and H’ are obtained from the block matrix defined in (4.11) 


To find A, 


(4.17) 


is enforced. 


If (4.12) is substituted into (4.17), 


Nac aC 
cae llr=ie) ae PENG, j= 75) late 


is obtained. 


sail 


dc , 
ae | et and =A 
r=1],6 


or. oor 





Allow - a vet ij 


(4.18) 


(4.19) 


Finally from (4.18) and (4.19), H, =D ; occurs and we solve for A; as A = AvH. 


These A values can now be used in Equation (4.12) to arrive at the values of w and c. 


V. RESULTS AND DISCUSSION 


Numerical calculations were performed for following parameters. 
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Table 1. Parameters Used in Numerical Calculations 


One of the main problems associated with many investigations that deal with 
complete numerical solution of Navier Stokes equations is the behavior of numerical 
scheme at high Reynolds numbers. It was found here that numerical results could be 
obtained for higher values of Reynolds numbers. For computational purposes, the results 


for Reynolds numbers up to 5000 are presented. 


For comparison, T (Torque) is defined as 
T= |o',r’Sindds’ (5.1) 
S,r=) 


If (5.1) is nondimensionalized, the following occurs: 


T 


2na°pvV C2 


r 
=[ogSit0d0 whee 0 «12 (te) ade Ye 
0 

The smal! mapping parameter, L, with increasing Reynolds number provides a finer 
grid in the boundary layer. The re-circulation of the flow near the equator is another 
important result of a spinning sphere which will be explained later. Since most of the sign 
changes in velocity profiles occur in this region, decreasing L much below one is not a 
solution to gain some additional grid points in the boundary layer for high values of 
Reynolds number because finer grid near the sphere means coarser grid in the outer 
region. 

Several solutions for large Reynolds numbers have been obtained by integrating 


numerically the boundary layer equations found in [Ref. 13] which gives 


6.48 a°@ 
Nis WNC hc agit: — lames (5335 
Re |” * Vv 


There are some other investigations that recommend a different numerical factor 
such as 5.95 [Ref. 8], 6.54 [Ref. 12] and 6.53 [Ref. 14]. In [Ref. 26], Equation (5.3) was 


modified as 


6.48 3] 


—— + 5.4 
Re,” Re, On 


The Constants 6.48 and 31 in Equation (5.4) are for Reynolds numbers based on radius; 
these constants should be modified to 9.16 and 62 respectively for Reynolds numbers 


based on diameter as in this study. 
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Table 2. Nondimensional Torque (M) at Different Reynolds Numbers 















Present | 








Figure | shows the variation of nondimensional torque with Reynolds numbers. It 
is seen that as the Reynolds number increases the general trend is quite consistent with the 


boundary layer solution given in [Ref. 13]. 





Figure 1. ---, Banks [Ref. 13]; 0, Sawatzki [Ref. 19]; +, Present 
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Figure 2 and Figure 3 show the transverse and azimuthal components of the 
dimensionless skin friction on the sphere, respectively. As it is in Figure 1, the higher 
Reynolds number, the closer it 1s to the boundary layer solution. Transverse and azimuthal 


components of dimensionless skin friction are defined as 


-(= ) Ug 
T= Re >. _ (5.5) 


7 2 (2) ; 
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Figure 2. Transverse Component of Skin Friction ; ---, Banks [Ref. 13] 
0, Rey=20; +, Rey=200; *, Reg=3000 
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Figure 3. Azimuthal Component of Skin Friction; ---, Banks [Ref. 13] 
0, Reg=20; +, Reg=200; *, Reg=3000 

Present results agree well with the boundary layer solution near the pole but this is 
not so for the equator region. This is because the boundary layer solution is obtained by 
integrating parabolic partial differential equations from 8=0 with known initial conditions 
and the solution near the equator is determined from this procedure. Thus, it does not 
satisfy the physical boundary conditions. 

Since a finite r..is used and the modified stream function is forced to be zero at r.. 
and on the sphere, a spurious re-circulating region of extremely weak flow (relative 
magnitude = 10”) is introduced to account for the continuity in the flow domain. The 
inflow at the pole changes to an outflow at the equator. The angular position of this 
change depends on Re but does not vary greatly in terms of radial distance from the 
sphere. However, the general trend is to get closer to the sphere with increasing Reynolds 
numbers. Figure 4 shows the stream lines in symmetric plane for different values of 


Reynolds numbers. 


ae) 
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Re=3000 
Figure 4. Stream Lines for Various Re 


In Figure 5, the jet effect can be easily seen at the equator for two different 
Reynolds numbers. The darker color represents the regions where the magnitude of 
dimensionless velocity in the symmetric plane is higher. The reference velocity in Figure 5 


is the velocity in the azimuthal direction on the sphere which has a value of unity. 
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Figure 5. Jet Effect at Equator 
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Figure 6a and Figure 6b show the experimental flow visualization obtained by 
F. P. Bowden and R. G. Lord [Ref. 17], and present numerical solution for Reg=3300 
respectively. Note that, unlike Figure5, the velocity fields in Figure 6a and Figure 6b are in 
3-D. As it is seen there is a reasonable resemblance between experimental and numerically 
obtained pictures. In the experiment, the sphere was levitated by means of a magnetic 
field. It was noted in [Ref. 17] that, because of this magnetic field, there is a slight upward 
convection current so the smoke rises slowly; this explains why the top half of the picture 
is white and the bottom half black. Note also that, the length of the radial jet is relatively 
shorter in Figure 6a. This behavior can be explained by the presence of possible shear 
turbulence in the radial jet for high Reynolds numbers causing the laminar jet behavior to 
disintegrate. Since the present results are obtained by means of a numerical solution 
without any kind of turbulence modeling, there is no effect of turbulence in our cases, 


even for higher Reynolds numbers. 





Figure 6a. Smoke Picture of the Boundary Layer Flow for Re=3300. (From Ref. 17) 





Figure 6b. Numerical Picture of the Boundary Layer Flow for Re=3300. 
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In the boundary layer approximation, it 1s assumed that the velocity profile in the 
boundary layer just before the equator 1s the same as the one at the start of boundary 
collision. This assumption yields zero radial velocity at 89 = 7/2, near the sphere. From 
the results it is seen that radial velocity is not zero at the equator near the sphere. Figure 7 
shows the variation of radial velocity with 9 for Re = 200 at different radial distances from 
the sphere. Note that the maximum value of radial velocity occurs exactly at the equator. 
Besides, there 1s an inward flow at both poles. Although the maximum value of radial 
velocity at the equator increases with Reynolds numbers, the thickness of the region where 


a general outward flow from the equator can be discussed decreases. 
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Figure 7. Variation of Radial Velocity with 9 for Re = 200 
-.-., f=] .0543; ---, r=1.1129;__, r=1.1651 


Figure 8 shows the variation of radial velocity with radial distance at the equator 
for different values of Reynolds numbers. Details of this swirling jet at the equator were 


given in [Ref. 27] and [Ref. 28]. 
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Figure 8. Variation of Radial Velocity with Radial Distance at @=1/2 0, Re=20 
+, Re=200; *, Re=3000 
Figure 9 shows the variation of dimensionless vorticity at the surface with © at 
different Reynolds numbers. The interest of this figure is to show that no separation 
occurs over the surface of the sphere, since a sign change must take place as a point of 


separation is passed. 





A) 


Figure 9. Variation of Dimensionless Vorticity with 0 
0, Re=20; +, Re=200; *, Re=3000 


Zi 


It 1s known that the boundary layer solution to be approached as Reynolds 
numbers increases, especially near the poles, and it is accurate to O(Re’~). Therefore, we 


da 
might expect that Re" 59 ars | and 6 = 0 will be of the form 


] da) ayn 
eee : (5.7) 


where A and B are constants. 


OW 
If Re~’? — at r=1 and 0 = 77/2 is evaluted, it can be seen that this function is an 


increasing function of Re. It can be assumed that this expression is in the form of 


suacae =CRe ° 5.8 
00 sae - . | 





where C and @ are constants. 


The values of A, B, C and @ were found by Dennis, Ingham and Singh [Ref. 26] to 


be 0.51, 0.56, 0.6 and % respectively. The modified values of B and C for Reynolds 
numbers based on diameter are 0.79 and 0.5 respectively 





, OW 
Table 3. Calculated Values of Re~”’? — atr=1,0=QOandr=1,0=7/2 
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VI. CONCLUSION AND RECOMMENDATIONS 


In addition for the hmiting cases of Re<<! and Re>>1 for which Stokes flow and 
boundary layer solutions are available, the numerical results showed very good agreement 
with the those previously obtained by other numerical techniques. By changing some of 
the parameters such as r.. and L, the agreement with previous numerical results could be 
refined, although there is no proof that the other studies are correct within three decimal 
digits. | 

Choosing the mapping parameter, L, should be done carefully so as to have many 
grid points in the regions of large gradients. A much more powerful mapping should be 
determined that provides a finer grid in two completely different regions, namely the 
boundary layer region and the re-circulation region (which though spurious, must be 
captured accurately by the numerical method for other calculations). 

A different FFT technique can be used that allows a number of grid points in the 8 
direction other than the ones that are a power of two. This would provide a relaxation in 
the choice of the number of grid points for different values of the Reynolds number which 
greatly affect the memory requirement. The problem is solved for 8 = 0 to 8 = & with the 
aid of sine expansions. Since the problem is also symmetric with respect to 7/2, by taking 
care of the even function properties it can only be solved for one quarter which simply 
doubles the grid points in 9 direction with the same memory requirement. 

The Green’s Function Method that deals with the Neumann boundary conditions 
should not be considered a trivial detail and must be used to have a stable convergence 
process. 

Obtained modified stream function values, C, can be easily used to solve the heat 
transfer problem due to convection from a spinning sphere. Derived additional energy 


equations are as follow 


De) 


27)? Laat —— 21r2 re an - 
reat) D Be) re) °D + Bo + Be [SKU K(t-Ar)] (6.1) 


where Ec = Pr Re. 
The nonlinear term, K, in Equation (6.1) consists of 5 parts 


_ A AC _o oF oC 
36) Bae 26 ae 
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Since, for forced convection only, Equation (6.1) is uncoupled from Navier Stokes 
equations, the steady state values of modified stream function, C, can be used in the 
energy equation after solving the fluid problem which decreases the additional memory 


requirement due to different block matrixes in Equation (6.1). 
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APPENDIX A. PROGRAM STRUCTURE 


In the first part of the program, problem parameters were defined such as Re, dt, 
L, r_inf, M, N. Also calculations that are not dependent on time were performed in this 
part. To gain some additional memory D1, D2 and D3 were calculated with D_Bmat 
function and saved as diagonal elements only. Non-diagonal elements of these matrixes 
were included in each iteration. 

Calculations showed that the matrix M1 (derivative operator in r direction) and 
consequently M3 and block matrixes that are diagonals of D1, D2 and D3 are ill- 
conditioned. To avoid numerical errors in solving equation systems, the Singular Value 
Decomposition method was used. Necessary SVD decomposition of the matrixes and 
calculations of w and c that are used in the Green’s Function Method were performed in 
this part. The built-in function in MATLAB, svd, was used for decomposition. 

Time dependent calculations basically consist of four steps. In the first step, Ue 
was calculated from Equation (3.32). In the second step, particular solutions of @ and c 
were found from Equations (3.25) and (3.26). In the third step, H, which is used in the 
Green’s Function Method, was found and finally homogenous and the particular solutions 
were combined to get wand c. 

Most of the calculations were performed in spectral domain which is a feature of 
the numerical technique. Only nonlinear operations were done in_ physical domain. 
Nonlinear terms were calculated with the functions nonlin and nonlinsp for @c and u 
equations respectively. 

A convergence criterion was not used in the program. For high Reynolds Number 
(above 200), some over-shooting occurs before steady state. This situation might have 
caused wrong results if a criterion was used. Instead, the results were checked during 


calculations and the program was terminated after seeing steady state. 


oF 


Solsvd : This function was used to solve equations after SVD decomposition. The 
limit used in this function reduces numerical errors due to ill-condition matrixes by 
skipping close to singular values in each iteration. 

Solsvd] : Since this function was used for solving homogenous parts in the 
Green’s Function Method and performed once before the time-loop, no limit was used in 
solsvd]. 

Phy : This function was used for transformation of the variables from spectral 
domain to physical domain 

Operat : Derivatives in 8 direction in the spectral domain were performed with the 
function operat. Since these terms were used only in nonlinear terms, the output from this 
function is in the physical domain. 

Variables that were used in the program are 

L : Mapping parameter. 

Re — : Reynolds Number. 

dt: Time interval between iterations. 

r_inf : Outer boundary. 

step : Number of iterations. 

N —_: Number of grid points in 6 direction. 

M — : Number of grid points in r direction. 

xcor, ycor : Cartesian coordinates of grid points with respect to the center of the 

sphere (used for the presentation of the results). 

M1 _ : Derivative operator in r direction. 


D1, D2, D3 : Matrixes defined as in Equation (4.11). 


Ww : Vorticity. 
c : Potential function as in Equation (3.15). 
u : Velocity component in ® direction. 


wtil, ctil : Homogenous part of the @ and c defined as in the Green’s Function 
Method. 


Aij, H, landa : Variables defined as in the Green’s Function Method. 
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wp, cp : Particular solution of @ and c defined as in the Green’s Function Method. 


torq  : Dimensionless torque defined as in Equation (5.2). 


5), 





APPENDIX B. PROGRAM CODES 


%% define program parameters 


clear all] 

global E M1 M2 M3 
L=7:Re=2-4t=0,05: 
r_inf=150;step=4500; 


%% define grid points 


N=32;n=[0:N-1]'; 
teta=n.*pi./N; 
M=32;m=[0:M]’; 
Z=cos(p1.*m./M); 
b=1+2*L/(r_inf-1); 
r=1+L.*(1+z)./(b-z); 
global NM 
xcor=cos(teta)*r’; 
ycor=sin(teta)*r’; 


E=emat(M); 
disp(E matrix is done’) 


%% define M1,M2,M3 


dz_dr=(b.*(L+r-1)-b.*(r-1)+L)./(L+r-1).42; 
B=diag(dz_dr); 

M1=B*E; 

R=diag(r);global R Re dt 

M2=R*M1; 
M3=R*(2.*eye(M+1)+M2)*M1; 


%% define D1 D2 D3 
(D1,D2,D3J=D_Bmat(M,N); 

% define initial wc wpcpu 
wp=zeros(M+1,N);cp=wp;w=wp;c=cp; 


dum 1=sin(teta)*(1./r.*2)';dum2=[dum 1 ;zeros(1,M+1);-1.*flipud(dum1(2:N,:))]; 
dum3=fft(dum2); 
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u=wp; 


% define kmat (for teta der.) and cotteta (for nonlinear terms) 


dum=ones(1,M+1); 
kmat=n*dum; 
cotteta=[0;cot(teta(2:N))]*dum; 
global kmat cotteta teta 


% svd decomposition for D1 D3 and define rsquare (for nonlinear terms) 


D1lu=zeros(N*(M+1),M+1); 
D11=D1u;D1p=D1u;D3l=D1u;D3u=D1u;D3p=D Lu; 
rsquare=zeros(M+1,N); 
for 1=1:N 
bas=(i-1)*(M+1)+1; 
son=bas+M;; 
rsquare(:,1)=1.“2; 
D1(bas,:)=[1 zeros(1,M)]; % modify D1 for w @ r=inf 
D1(son,:)=[zeros(1,M) 1]; % modify D1 for w @ r=] 
{D11(bas:son,:),D1u(bas:son,:),D 1 p(bas:son,:)}=svd(D 1 (bas:son,:)); 
D3(bas,:)=[1 zeros(1,M)]; % modify D3 for c @ r=inf 
D3(son,:)=[zeros(1,M) 1]; % modify D3 forc @ r=1 
[D31(bas:son,:),D3u(bas:son,:),D3p(bas:son,:) ]=svd(D3(bas:son,:)); 
end 
global rsquare 


% define c_c_tetac_ru_u_tetau_r 
Ce=Zeros( INV |) Cateta=c. Car=c_: 
Hele simetetay. (l/r 2 i: 

Weteta—1. cos(tetay (l/r. 2): 
Met——o a sHmtetal (lasiee 5)! 
Slopalice catera can leu tela War 

% define c boundary at r=1 
cbr1=zeros(1,N); 


% begin green function wtil ctil 


wtil=zeros(M+1,N“%2);ctil=wtil; 
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for jJj=1:N 
RR=zeros(M+1,N); 
dum1=[zeros(1,jj-1) 1 zeros(1,N-jj)]; 
dum2=[dum]1 0 -1.*flipir(dum1(2:N))]; 
dum3=fft(dum2); 
dum4=imag(dum3(1:N)); 
RR(M+1,:)=dum4; 
wtl(:,jj*N)=solsvd1(D 11((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:).... 
D1lu((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,,:).... 
D1 p((N-1)*(M+1)+1:(N-1)*(M+1)+1+M.,:).... 
RR(:,N)); 


for 1=N-1:-1:1 
sum=zeros(M+1,1); 
11=(jj-1)*N-+H; 
for j=i+1:N 
if rem(i+j,2)==2 
sum=sum+(-2*(i-1)).*wtil(:,Qj-1)*N-+)); 
end 
sum(1)=0; 
sum(M+1)=0; 
end 
RR_=RR(-:,1)-sum; 


wtil(:,11)=solsvd1(D11((1-1)*(M+1)+1:G-1)*(M+1)+1+M,;:).... 
D 1u((1-1)*(M+1)+1:(1-1)*(M+1)+14M,:).... 
Di p((i-1)*(M+1)+1:(C-1)*(M+1)+1+M,:).... 
RR_); 


end 


tempwtil=(-1.*rsquare).*wtil(:,Qj-1)*N+1:jj*N); 

tempwtil(M+1,:)=zeros(1,N); 

tempwti(1,:)=zeros(1,N); 

ctil(:,jj*N)=solsvd1(D31((N-1)*(M+1)+1:(N-1)*(M+1)+1+M.:).... 
D3u((N-1)*(M+1)+1:(N-1)*(M+1)+1+M.,:).... 
D3p((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:).... 
tempwtil(:,N)); 


for 1=N-1:-1:1 
li=(jj-1)*N+; 
sum=zeros(M+1,1); 
for j=i+1:N 
if rem(i+j,2)==2 
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sum=sum-+(-2*(i-1)).*ctil(:,Qj-1)*N-+)); 
end 
sum(1)=0; 
sum(M+1)=0; 
end 
wtil_=tempwtil(:,1)-sum; 


ctil(:,1)=solsvd1(D31((1-1)*(M+1)+1:q-1)*(M+1)+1+M.,:),... 
D3u((i-1)*(M+1)+1:(-1)*(M+1)+14+M.,:).... 
D3p((i-1)*(M+1)+1:G-1)*(M+1)+1+M,:).... 
wtil_); 
end 


end 
% find Aij matrix 


for Jj=1:N 
blok=ctil(:,qj-1)*N+1:j*N); 
dum 1!=phy(blok’)’; 
dumlZ—=(Vile dumb); 
dum3=[dum2; zeros(1,M+1) ; -].*flipud(dum2(2:N,:))]; 
dum4=fft(dum3);dum5=imag(dum4(1:N,:));dum5=dum5'; 
ctilr_1(jj, 1:N)=dum5(M-+1,:); 

end 


Ay=ctilr_1’; 
% find svd of Ajj for landa solving 


Aij(1,:)=[1 zeros(1,N-1)]; 
[Aiju,Aijs, Aijv]=svd(Ajj); 


% begin time loop 
for ops=1:step 
if ops==1;prc=c;prw=w;pru=u;end 
% find R1 
Rl=zeros(M+1,N); 
R1¢,N)=D2((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:)*w(,N); 


for i=N-1:-1:1 
summ=zeros(M+1,1); 
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for j=i+1:N 
if rem(i+j,2)==2 
summ=summ-+(-2*(i-1)).*w(,J); 
end 
end 
R1(¢:,)=D2((i-1)*(M+1)+1:€1-1)*(M+1)+14+M,:)*w(.,1)+summ; 
end 
R1l=-1.*R1; 


RIisp=zeros(M+1,N); 
Ri sp(:,N)=D2((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:)*u(:,N); 
for 1=N-1:-1:1 
summ=zeros(M+1,1); 
for j=1+1:N 
if rem(i+j,2)==2 
summ=summ+(-2*(1-1)).*u(:,J); 
end 
end 
R1sp(:,1)=D2((i- 1)*(M+1)+1:(-1)*(M+1)+1+M,:)*uC,1)+summ; 
end 
Rlisp=-1.*R1sp; 


% find R2sp 
R2asp=nonlinsp(u,c); 
R2bsp=nonlinsp(pru,prc); 
R2sp=3/2.*R2asp-1/2.*R2bsp; 
R2sp=(-1*Re).*rsquare.*R2sp; 
Rtotsp=RIsp+R2sp; 

% find new u 


%% apply BC for u 


Rtotsp(1,:)=zeros(1,N); 9% BC for u=0 at r_inf 
Rtotsp(M+1,:)=zeros(1,N); % BC for u at r=1 


u(:,N)=solsvd(D11((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:).... 
D1lu((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:),... 
D1 p((N-1)*(M+1)+1:(N-1)*(M+1)+14M,;:),... 
Rtotsp(:,N)); 


for i=N-1:-1:1 


by, 


sum=zeros(M+1,1); 

for j=1+1:N 
if rem(i+j,2)==2 

sum=sum+(-2*(1-1)).*uC:,J); 

end 

sum(])=0; 

sum(M-+1 )=0; 

end 

Rtotsp_=Rtotsp(:,1)-sum; 


u(:,1)=solsvd(D 11((1-1)*(M-+1)+1:(i- 1)*(M+1)+14+M,:),... 
D1u(di-1)*(M+1)+1:-1)*(M+1)+14+M,:),... 
D1 p(di-1)*(M+1)+1:q-1)*(M+1)+14+M,:),... 
Rtotsp_); 


end 
% find R2 


R2a=nonlin(w,c,u); 
R2b=nonlin(prw,prc,pru); 
Re=3/2e 4-1/2, 2b, 
Re—(-lahe .rsquane, i 2: 
Rtot=R1+R2; 
prw=w;prc=c;pru=u; 


%% green function solution 
% find new wp 
%% apply BC for wp 


Rtot(1,:)=zeros(1,N); % BC for wp=0 at r_inf 
Rtot(M+1,:)=zeros(1,N); % BC for wp=0 at r_] 


wp(:,N)=solsvd(D11((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:).... 
D1lu((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:),... 
D1 p((N-1)*(M+1)+1:(N-1)*(M+1)4+14+M,:),... 
Rtot(:,N)); 


for 1=N-1:-1:1 
sum=zeros(M+1,1); 
for j=1+1:N 
if rem(i+j,2)==2 
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sum=sum+(-2*(i-1)).*wp(:,j); 
end 
sum(1)=0; 
sum(M+1)=0; 
end 
Rtot_=Rtot(:,1)-sum; 


wp(:,1)=solsvd(D 11((i- 1 )*(M+1)+1:(1-1)*(M+1)+14+M,:),... 
D1lu((i-1)*(M+1)+1:/-1)*(M+1)+14+M,:).... 
D1 p((i- 1)*(M+1)+1:(i-1)*(M+1)+14M,;),... 
Rtot_); 


end 
% solve for new cp 


tempwp=-1.*rsquare.* wp; 
tempwp(M-+1,:)=cbr1; % BC for cp=-c_r at r_] zero for pure spin 
tempwp(1,:)=zeros(1,N); % BC for cp=0 at r_inf 


cp(:,N)=solsvd(D31((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:).... 
D3u((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:),... 
D3p((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:).... 
tempwp(:,N)); 


for 1I=N-1:-1:1 

sum=zeros(M+1,1); 

for j=i+1:N 
if rem(i+j,2)==2 

sum=sum+(-2*(i-1)).*cp(:,j); 

end 

sum(1)=0; 

sum(M-+1)=0; 

end 

w_=tempwp(:,1)-sum; 


cp(:,)=solsvd(D31((i-1)*(M+1)+1:G-1)*(M+1)+14M.:).... 
D3u((i- 1)*(M+1)+1:(i-1)*(M+1)+14+M,°).... 
D3p((i-1)*(M+1)+1:(7-1)*(M+1)+1+M,:).... 
W_); 
end 


Jo find landa 
% find H 
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dum 1=phy(cp’)’; 

dum2=(M1*dum1)'; % NxM+1 

dum3=-1.*c_r-dum?; 
dum4=[dum3;zeros(1,M+1);-1.*flipud(dum3(2:N,:))]; 
dum5=fft(dum4); 

dum6=imag(dum5(1:N,:)); 

dum6=dum6'; % M+1xN 


H=dum6(M+1,:);H=H'; % Nx! 
ath uP 
landa=solsvd(Aiju,Atjs,Aljv,H); 


% find complete w & c 


dum 1!c=zeros(M+1,N);dum 1 w=dum Ic; 

for Jj=1:N 
blokc=ctil(:,Qj-1)*N+1:j*N); 
blokw=wtil(:,(jj-1)*N+1:jj*N); 
dum2c=landa(jj).*blokc; 
dum2w=landa(jj).*blokw; 
dum1]c=dum 1c+dum2c; 
dum 1w=dum1w+dum2w;: 

end 


w=wp+dum | w; 
c=cp+dum|c; 


dum | =(phy(u')+u_)'; 

dum2=M 1] *dum1-dum1; 
dum3=dum2(M-+1,:)'.*(sin(teta).*2); 
torq(ops)=-8*pi/Re*trapz(teta,dum3); 
if ops>1 

err(ops)=torq(ops)-torq(ops- 1); 

end 


save tez4.mat uc w Re ops r_inf torq LN M dt; 
end % end of time loop 


function [D1,D2,D3]=D_Bmat(M,N) 
global M3 R dt Re 

k=(M-+1)*N; 

D1=zeros(k,M+1); 
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D2=zeros(k,M+1); 
D3=zeros(k,M+1); 
for 1i=1:N 
fon j= 
il 
D1(C1-1)*(M+1)+1:(i-1)*(M+1)+1+M,:)... 
=M3+(-1*(-1)*(1).*eye(M+1)-Re/dt.*R.%2); 
D2((i-1)*(M+1)+1:(C/-1)*(M+1)+14+M.,,:)... 
=M3+(-1*G-1)*(1).*eye(M+1)+Re/dt.*R.%2); 


D3((i-1)*(M-+1 )+1:G-1)*(M+1)+14M,:)... 
=M3+(-1*(i-1)*(1).*eye(M+1)); 
end 
end 
end 


function E=emat(M) 
E=zeros(M,M); 
m=0:M; 
z=cos(pi.*m./M); 
for 1=0:M 
for j=0:M 
if i~=] 
if ((i==0 | i==M) & G==01 j==M)) | (i~=0 & i~=M) & (j~=0 & j~=M)) 
c=] ;else 
if G==0 | 1==M) & (G~=0 & j~=M) 
c=2;else 
if (i~=0 & 1~=M) & G==0 | j==M) 
c=1/2;end;end;end; 
E(i+1 jt 1)=c*(-1)4C4+))/(zi+ 1)-zg+1)); 
else; 
ifa== 
E(i+1j+1)=(2*M%2+1)/6; 
else 
ini== 
E(i+1,j+1)=(2*M‘%2+1)/(-6); 
else 
E(i+ 1 ,j+1)=-1*2G+1)/(2*(1-z9G+1)%2)); 
end 
end 
end 
end 
end 
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function R2=nonlin(w,c,u) 

global kmat rsquare cotteta NM M1 c_c_teta c_r teta u_ u_r u_teta 
% variable list 

Jo 

% dw/dteta=non! dc/dteta=non2 dw/dr=non3 dc/dr=non4 (in phy) 
% du/dteta=non5 du/dr=non6 

% transpose matrixes for colomn operations in teta direction 
r=sqrt(rsquare)';w=w';c=c';u=u'; 

non!l=operat(w); 

non2=operat(c); 

non5=operat(u); 

wph=phy(w);cph=phy(c);uph=phy(u); 

non3=[{[M1*wph'J; 

non4=[M 1*cph']';non4(:,M+1)=-0.5.*sin(teta); 

non6=[M 1 *uph']’; 

term1|=-1./r.*non3.*(non2+c_teta); 
term2=1./r.*non1!.*(non4+c_r); 
term3=-1.*cotteta./r.*wph.*(non4+c_r); 
term4=-1!.*cotteta./r.*non3.*(c_+cph); 

term5=wph./rsquare .*(non2+c_teta); 
term6=(c_+cph)./rsquare’.*non 1; 
term7=1.*cotteta.*(non6+u_r).*(uph+u_).*2./r; 
term8=-1.*(non5+u_teta).*2.*(uph+u_)./rsquare’: 

term0=term | +term2+term3+term4+term5+term6+term7+term8s; 
ter=[term0;zeros(1,M+1);-1.*flipud(term0(2:N,:))]; 

RR 2=fititer): 

R2=imaegeRR2(1:N,:));R2=R2’; 


function R2sp=nonlinsp(u,c) 

global kmat rsquare cotteta N M M1 c_c_teta c_r teta u_ u_r u_teta 
% variable list 

%o 

% du/dteta=non1 dc/dteta=non2 du/dr=non3 dc/dr=non4 (in phy) 
% transpose matrixes for column operations in teta direction 
r=sqrt(rsquare) ;u=u';c=C’; 

non!=operat(u); 

non2=operat(c); 

uph=phy(u);cph=phy(c); 

non3=[M1*uph'J; 

non4=[M 1*cph']’; 

term 1=-1./r.*(non3+u_r).*(non2+c_teta); 
term2=]1./r.*(non|+u_teta).*(non4+c_r); 
term3=cotteta./r.*(uph+u_).*(non4+c_r); 

term4=-] .*cotteta./r.*(non3+u_r).*(c_+cph); 
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term5=-1.*(uph+u_)./rsquare’.*(non2+c_teta); 
term6=(c_+cph)./rsquare'.*(non1+u_teta); 
term0=term 1 +term2+term3+term4+term5+term6; 
ter=[term0;zeros(1,M+1);-1.*flipud(term0(2:N,:))]; 
RRZ=fit(ter): 

R2=imag(RR2(1:N,:));R2sp=R2'; 


function y=operat(x) 

global N M kmat 

re=zeros(N,M+1); 

X=fese xX: 

=e imac 
Y_ful=[Y;zeros(1,M+1);flipud(Y(2:N,:))]; 
y=ifft(Y_ful); 

y=real(y(1:N,:)); 


function y=phy(x) 

global NM 

re=zeros(size(x)); 

X=ret+i.*x; 
Y_ful=[X;zeros(1,M+1);conj(flipud(X(2:N,:)))]; 
dum=ifft(Y_ful); 

y=real(dum(1:N,:)); 


function y=solsvd(u,s,v,b) 
w=diag(s); 
wmin=max(w)* le-12; 
for 1=1:length(w) 
if w(i)<wmin;ww(1)=0;else;ww(1)=1/w(i);end 
end 
www=diag(ww); 
y=v*www*(u'*b); 


function y=solsvd1(u,s,v,b) 

w=diag(s); 

for 1=1:length(w) 
ww(i)=1/w(i); 

end 

www=diag(ww); 
=v*www*(u'*b); 
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